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A Eulerian-Lagrangian  model  to  simulate 
two-phase/particulate  flows 

By  S.  V.  Apte,  K.  Maheshj,  & T.  Lundgren^ 


1.  Motivation  and  Objectives 

Figure  1 shows  a snapshot  of  liquid  fuel  spray  coming  out  of  an  injector  nozzle  in 
a realistic  gas-turbine  combustor.  Here  the  spray  atomization  was  simulated  using  a 
stochastic  secondary  breakup  model  (Apte  et  al.  2003a)  with  point-particle  approxima- 
tion for  the  droplets.  Very  close  to  the  injector,  it  is  observed  that  the  spray  density  is 
large  and  the  droplets  cannot  be  treated  as  point-particles.  The  volume  displaced  by  the 
liquid  in  this  region  is  significant  and  can  alter  the  gas-phase  flow  and  spray  evolution. 
In  order  to  address  this  issue,  one  can  compute  the  dense  spray  regime  by  an  Eulerian- 
Eulerian  technique  using  advanced  interface  tracking/level-set  methods  (Sussman  et  al. 
1994;  Tryggvason  et  al.  2001;  Herrmann  2003).  This,  however,  is  computationally  in- 
tensive and  may  not  be  viable  in  realistic  complex  configurations.  We  therefore  plan  to 
develop  a methodology  based  on  Eulerian-Lagrangian  technique  which  will  allow  us  to 
capture  the  essential  features  of  primary  atomization  using  models  to  capture  interac- 
tions between  the  fluid  and  droplets  and  which  can  be  directly  applied  to  the  standard 
atomization  models  used  in  practice.  The  numerical  scheme  for  unstructured  grids  de- 
veloped by  Mahesh  et  al.  (2003)  for  incompressible  flows  is  modified  to  take  into  account 
the  droplet  volume  fraction.  The  numerical  framework  is  directly  applicable  to  realistic 
combustor  geometries. 

Our  main  objectives  in  this  work  are: 

• Develop  a numerical  formulation  based  on  Eulerian-Lagrangian  techniques  with 
models  for  interaction  terms  between  the  fluid  and  particles  to  capture  the  Kelvin- 
Helmholtz  type  instabilities  observed  during  primary  atomization. 

• Validate  this  technique  for  various  two-phase  and  particulate  flows. 

• Assess  its  applicability  to  capture  primary  atomization  of  liquid  jets  in  conj  action 
with  secondary  atomization  models. 


2.  Mathematical  Formulation 

Recent  direct  numerical  simulations  of  large  number  of  solid  particles  interacting 
through  a fluid  medium  by  Joseph  and  collaborators  (Choi  & Joseph  2001)  show  that  a 
layer  of  heavy  particles  with  fluid  streaming  above  it  can  develop  Kelvin-Helmholtz  (K- 
H)  instability  waves  whereas  a layer  of  particles  above  a lighter  fluid  develops  Rayleigh- 
Taylor  instability.  This  suggests  that  the  primary  breakup  of  a liquid  jet  into  a spray  can 
be  simulated  by  replacing  the  jet  by  a closely  packed  collection  of  droplets  with  some 
assumed  size  distribution.  The  K-H  instability  at  the  boundary  between  droplets  and 
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Figure  1 . Snapshot  spray  from  a gas-turbine  fuel-injector. 

fluid  would  initiate  dispersal  of  droplets  into  a spray.  Further  breakup  of  these  dispersed 
drops  can  be  obtained  by  advanced  secondary  breakup  models. 

The  formulation  described  below  is  a modification  of  the  equations  for  spray  compu- 
tation developed  by  Dukowicz  (1980)  which  consists  of  Eulerian  fluid  and  Lagrangian 
particle  calculations,  and  accounts  for  the  displacement  of  the  fluid  by  the  particles  as 
well  as  the  momentum  interchange  between  them.  The  modifications  presented  here  are 
mainly  in  the  details  of  modeling  the  interaction  terms. 

2.1.  Gas-Phase  Equations 

The  fluid  mass  for  unit  volume  satisfies  a continuity  equation, 

Q 

^ (p/0/)  + V • = 0 (2.1) 

where  p/,  0/,  and  u/  are  fluid  density,  volume  fraction,  and  velocity,  respectively.  This 
indicates  that  the  average  velocity  field  of  the  fluid  phase  does  not  satisfy  the  divergence- 
free  condition  even  if  we  consider  an  incompressible  suspending  fluid.  The  fluid  momen- 
tum equation  is  given  as 

^ (P/0/U/)  -b  V • (/’/0/U/U/)  = - V (0/P)  + V • (/^/Dc)  + F (2.2) 

where  p is  the  average  dynamic  pressure  in  the  fluid  phase,  p/  is  the  viscosity  of  the 
fluid  phase,  and  Dc  = V^c  + is  the  average  deformation-rate  of  the  fluid-particle 
composite,  Uc  is  the  composite  velocity  of  the  mixture,  and  F is  the  force  per  unit 
volume  the  particles  exert  on  gas.  These  equations  axe  derived  in  detail  for  constant 
density  flows  by  Joseph  & Lundgren  (1990).  For  particulate  flows  and  dilute  suspensions 
at  low  Reynolds  numbers,  the  fluid  viscosity  should  be  replaced  by  an  effective  viscosity 
p*  by  using  Thomas  (1965)  correlation, 

p*  =pf{\  + 2.50/  -b  10.050^  -b  0.00273e^®'®®^)  (2.3) 

2.2.  Particle-Phase  Equations 

The  evolution  of  particle-phase  is  governed  by  a Liouville  equation  for  the  particle  dis- 
tribution function  $(Xp,  Up,pp,  V^,t) 
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— + Vx  • (^Up)  + Vup  ■ (^Ap)  = 0,  (2.4) 

where  Xp  is  the  particle  position,  Up  particle  velocity,  pp  particle  density,  and  Vp  particle 
volume.  Ap  is  the  particle  acceleration  and  Fp  = rUpAp  the  total  force  acting  on  the 
particle  of  mass  nip  and  are  given  below.  Here,  Vx’  and  Vup’  are  the  divergence  operators 
with  respect  to  space  and  velocity,  respectively.  The  individual  particle  positions  and 
velocities  can  be  obtained  by  solving  the  Liouville  equation  in  Lagrangian  framework  for 
each  particle  p: 


1 \ 

dt  <*'> = ">■ 

(2.5) 

rUp—  (up)  = Fp 

(2.6) 

2.2.1.  Particle  Forces 

The  main  issue  is  to  model  the  force  on  a particle.  This  may  consist  of  the  standard 
hydrod5mamic  drag  force,  dynamic  pressure  gradient,  gradient  of  viscous  stress  in  the  fluid 
phase,  a generalized  buoyancy  force,  and  inter-particle  collision.  The  total  acceleration 
of  the  particle  Ap  is  given  as, 

Ap  = Dp  (u/  - Up)  - -i-  VFp  + ( 1 - — J g + Bp  -t-  Acp  (2.7) 

Pp  \ Pp/ 

Here  Bp  is  the  generalized  buoyancy  force  and  Acp  is  the  acceleration  due  to  inter- 
particle forces.  If  Pp  >>  pg  the  pressure  gradient,  viscous,  and  buoyancy  terms  are 
usually  negligible.  In  the  present  study,  the  generalized  buoyancy  force  is  also  neglected 
for  simplicity.  It  is  shown  later  that,  even  without  the  presence  of  this  buoyancy  force, 
one  can  obtain  lift  of  particles  in  a shear  flow.  The  drag  force  is  caused  by  the  motion  of 
a particle  through  the  gas.  The  standard  expression  for  Dp  is  used 


Dr,  = ^C, 


PS  I"/  ~ "P 


(2.8) 


8 Pp  ± Lp 

where  is  the  drag  coefficient  and  is  given  by  (Gidaspow  1994;  Andrews  & O’Rourke 
1996) 


Cd  = — (1  + afle^)  Oji’®,  for  Rcp  < 1000  (2.9) 

= O.440J^  ®,  for  Rep  > 1000  (2.10) 

where  Ca  is  the  drag  coefficient  for  spherical  particles,  Rp  = {SVp/An)^^^  is  the  particle 
radius.  The  particle  Reynolds  number  {Rcp)  is  given  as 

Pf 

There  is  an  indirect  collective  eflFect  in  this  drag  term:  when  there  is  a dense  collection 
of  particles  passing  through  the  fluid  interphase  momentum  exchange  term  in  equation 
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(2.2)  will  cause  Ug  to  approach  the  particle  velocity,  Up  thus  decreasing  the  drag  on  a 
particle,  a drafting  effect. 

The  probability  function  $ integrated  over  velocity  and  mass  gives  the  probable  number 
of  particles  per  unit  volume  at  x and  t in  the  interval  (up,  Up  4- dup) , (pp,  pp  + dpp) , 
(^)  Vp  + dVp).  The  particle  volume  fraction  (0p)  is  defined  fi'om  the  particle  distribution 
function  as, 


Qp=  j j j ^VpdVpdppdUp 


(2.12) 


Prom  continuity,  the  gas-phase  volume  fraction  is  obtained  as  0^  = 1— ©p.  The  interphase 
momentum  transfer  function  per  unit  volume  in  equation  (2.2)  is  given  as 


F = //  |WpPp[z)p(u/-Up)-lv 


Pp 


dVpdppdvLp. 


(2.13) 


2.2.2.  Collision  Force 

The  acceleration  of  particles  due  to  inter-particle  interactions  (Acp)  is  an  important 
term  in  dense  particulate  and  two-phase  flows.  For  dilute  and  lightly  loaded  configura- 
tions, the  particle  volume  fraction  (0p)  is  small  (<  10%),  the  inter-particle  collisions  are 
negligible  and  probability  of  particles  overlapping  each  other  is  low.  For  dense  partic- 
ulate flows,  however,  the  particle  volume  fraction  should  not  exceed  the  close-packing 
limit  (which  is  usually  around  0.6  for  three-dimensional  case).  In  the  Eulerian-Eulerian 
approach  for  two-phase  flows,  this  is  ensured  by  force  due  to  the  gradient  of  interparti- 
cle stress  in  the  averaged  momentum  equation  for  the  particle  phase  (Gidaspow  1986; 
Gidaspow  1994).  Same  model  was  used  in  Eulerian-Lagrangian  approach  by  Andrews 
O’Rourke  (1996),  Patankar  k Joseph  (2001b),  Snider  et  al.  (1998),  and  Snider  (2001). 
Accordingly,  the  expression  for  acceleration  of  particles  due  to  collision  (Acp)  is 


€)pPp 


VT 


(2.14) 


where  r is  the  interparticle  stress  that  provides  a pressure  type  force  that  prevents  packing 
of  particles  beyond  the  close-packing  limit.  Expression  for  r is  given  as  (from  Snider  et 
al.  (1998)) 


T = 


©cp  0p 


(2.15) 


where  has  units  of  pressure,  ©cp  is  the  particle  volume  fraction  at  close  packing  and  ^ is 
a constant.  The  values  for  Ps  and  /?  are  obtained  from  Snider  et  al.  (1998).  In  this  model, 
the  inter-particle  acceleration  due  to  particle  collision  is  assumed  to  be  independent 
of  its  size  and  velocity.  This  model  is  least  expensive  in  terms  of  computational  cost 
as  particle  binary  pairs  are  not  formed  and  the  collision  force  is  directly  obtained  by 
interpolation  from  equation  (2.14).  This  model,  however,  couldn’t  completely  prevent  the 
particle  volume  fraction  to  exceed  the  close  packing  limit  and  some  numerical  instabilities 
were  experienced  in  the  present  unstructured  code.  Further  research  on  this  model  will 
be  conducted  to  eliminate  this  problem.  An  alternative  but  expensive  collision  scheme 
based  on  the  distinct  element  method  (DEM)  of  Cundall  & Strack  (1979).  This  scheme 
can  be  readily  applied  to  parcel  techniques  used  in  Lagrangian  spray  simulations.  A 
parcel  represents  a number  {Np)  of  droplets/paxticles  of  same  size,  velocity,  and  other 
properties.  The  effective  radius  {Rpar)  of  a parcel  based  on  its  total  volume  is  then  given 
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by 


(2.16) 


In  this  method,  the  interaction  among  particles  and  between  wall  and  particles  is  taken 
into  account  separately  to  ensure  that  0p  < ©cp.  The  force  on  parcel  p due  to 

collision  with  parcel  j is  given  by 


pP-P 


pj  — H for  dpj  ^ {PpanP  ^ Pparyj  “b 

= - Vc(  Up  Uj)  • Upj^  llpj  for  dpj  < {RparyP  "b  Rpar>j  "b  Oi) 

^Pj  ~ {PpafiP  "b  Rpari  j "b  tt)  dpj 


fpP-P 

■ jp 


= 

PJ 


p-p 


where  dpj  is  the  distance  between  the  center  of  the  and  parcels,  Upj  is  the  unit 
vector  from  the  center  of  parcel  j to  that  of  parcel  p,  a is  the  force  range,  the  stiff- 
ness parameter,  and  rjc  the  damping  parameter.  Tsuji  et  al.  (1993)  used  the  following 
expressions  to  compute  the  damping  parameter 


Vc 


I rupkc 
a — —In  (c/tt) 


where  e is  the  coefficient  of  restitution.  Similarly,  the  parcel- wall  force  on  parcel 

p due  to  collision  with  wall  w is 

= 0 dpw  > {RparyP  + Oi) 

“ ^c(^p)  ' ^pv)  for  dpiD  i^RpaV)  P “f" 

^pw  — {Rpar^P  "b  Ck)  dpyj 


where  dp^  is  the  distance  between  the  parcel  and  the  wall,  and  np„,  is  the  unit  vector  from 
the  wall  to  the  center  of  the  parcel.  The  total  colhsion  force  is  obtained  by  looping  over 
all  particles  and  walls.  The  corresponding  particle  acceleration  is  obtained  by  dividing 
the  collision  force  by  parcel  mass  (rUp  = NpPpVp). 

3.  Numerical  Method 

The  collocated,  finite-volume  numerical  scheme  on  unstructured  grids  developed  by 
Mahesh  et  al.  (2003)  is  modified  to  take  into  account  the  gas-phase  volume  fraction.  In 
addition,  the  particle  centroids  are  tracked  using  the  Lagrangian  framework  developed  by 
Apte  et  al.  (2003b).  The  important  feature  of  the  numerical  scheme  is  the  computation 
of  ©y  and  ©p  on  the  unstructured  grids.  As  the  particles  in  the  simulations  performed 
do  not  move  out  of  the  domain  or  are  destroyed,  the  total  volume  occupied  by  the 
particles  remains  constant  from  mass-conservation.  As  a particle  crosses  a particular 
control  volume,  its  contribution  to  the  particle  volume  fraction  on  neighboring  cells  must 
provide  global  mass-conservation.  In  addition,  it  was  observed  that  the  Eulerian  volume 
fraction  field  should  also  be  smooth  in  order  to  avoid  numerical  instabilities  encountered 
due  to  changes  in  © as  the  particles  move  from  one  grid  cell  to  another.  A strategy  similar 
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Computational  domain,  0.2  x 0.6  x 0.0275m 
Fluid  density,  \.2hAkg/m^ 

Numer  of  Parcels,  2880 
Diameter  of  particles,  500^m 


Grid,  10  X 30  X 5 
Particle  Density,  2500kg/m^ 
Particles  per  parcel,  3375 
Initial  particle  concentration,  0.4 


Table  1 . Parameter  description  for  the  gravity-dominated  sedimentation  case. 


Figure  2.  Temporal  evolution  particle  distribution  during  gravity-dominated  sedimentation. 
Initially  the  particles  are  randomly  distributed  over  the  box  rectangular  box. 


to  the  two-way  coupling  methodology  by  Maxey  & Patel  (2001)  is  used  to  compute  the 
volume  fraction.  The  interphase  momentum  transport  terms  are  also  treated  in  a similar 
way. 

The  particle  equations  are  integrated  using  third-order  Runge-Kutta  schemes  for  ode- 
solvers.  At  each  Runge-Kutta  step,  the  particles  were  located  and  the  collision  force 
was  re-computed.  This  was  found  necessary  to  ensure  that  the  close-packing  limit  for 
particle-volume  fraction  is  not  exceeeded. 


4.  Results 

We  first  simulated  sedimentation  of  solid  particles  under  gravity  in  a rectangular  box. 
Details  of  this  case  are  given  in  Table  1.  The  initial  parcel  positions  are  generated  ran- 
domly in  the  top  half  of  the  box,  0.3  < j/  < 0.6.  These  particles  are  then  allowed  to  settle 
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Computational  domain,  0.2  x 0.6  x 0.0275m 
Gas  jet  velocity,  9m/ s 
Fluid  density,  1.2Mkg/m^ 

Numer  of  Parcels,  2880 
Diameter  of  particles,  500fim 


Grid,  10  X 30  X 5 
Jet  diameter,  0.04m 
Particle  Density,  2500A:g/m® 
Particles  per  parcel,  3375 
Initial  particle  concentration,  0.4 


Table  2.  Parameter  description  for  the  simulation  of  fluidization  by  a gas  jet. 


through  the  gas-medium  under  gravity.  The  collision  model  is  important  here  near  the 
bottom  wall  after  the  initial  group  of  parcels  hit  and  bounce  back  from  the  bottom  wall. 
This  prevents  particle  volume  fraction  to  exceed  the  close-pack  limit.  The  upper  mixture 
interface  between  the  particles  and  the  fluid  is  closely  approximated  hy  h = gt'^ /2  similar 
to  that  obtained  by  Snider  (2001).  The  particles  eventually  settle  down  with  close-packing 
near  the  bottom  wall.  Snider  (2001)  and  Patankar  & Joseph  (2001a)  did  sedimentation 
under  gravity,  in  an  inclined  pipe  with  large  number  of  parcels  and  used  the  particle 
stress  model  for  collision  force.  This  allowed  them  to  use  large  number  of  parcels  as  the 
computational  cost  for  this  collision  model  is  negligible.  We  plan  to  use  this  model  to  do 
similar  validation  test  cases  using  the  unstructed  LES  code.  One  difiiculty  at  the  time 
of  writing  was  that  this  model  didn’t  always  prevent  the  particle  volume  fraction  from 
exceeding  the  close-pack  Umit.  Further  investigations  on  this  matter  will  be  performed. 

4.1.  Gas-Solid  Fluidization 

We  consider  the  problem  of  fluidization  of  solid  particles  arranged  in  an  array  at  the 
bottom  of  a rectangular  box.  Fluidization  is  achieved  by  a jet  of  gas  issuing  from  the 
bottom  of  the  box.  The  flow  parameters  are  given  in  Table  2.  The  choice  of  collision 
parameters  is  based  on  the  recent  computations  by  Patankar  & Joseph  (2001b).  The 
particle  motion  is  mostly  dominated  by  the  hydrodynamic  drag  force  and  collision  model 
should  not  affect  the  overall  particle  motion.  The  collision  model,  however,  is  important 
in  governing  the  particle  behavior  near  the  walls  and  helps  prevent  the  volume  fraction 
from  exceeding  the  close-pack  limit. 

Figure  3 shows  the  position  of  parcels  at  different  times  during  bubbling  fluidization. 
Parcel  diameters  are  drawn  to  scale.  The  jet  issuing  from  the  bottom  wall  pushes  the 
particles  away  from  the  center  region  and  creates  a gas-bubble  in  the  center.  The  particles 
collide  with  each  other  and  the  wall  and  are  pushed  back  towards  the  central  jet  along  the 
bottom  wall.  They  are  then  entrained  by  the  jet  and  are  levitated.  This  eventually  divides 
the  central  bubble  and  two  bubbles  are  trapped.  The  particles  tend  to  move  upwards  and 
collide  with  the  upper  wall  and  remain  levitated  during  future  times.  The  computational 
results  are  in  good  agreement  with  the  simulations  of  Patankar  & Joseph  (2001b)  as 
well  as  in  qualitative  agreement  with  experiments  on  jet  fluidization.  Similar  results 
are  reported  using  Eulerian-Eulerian  approach  in  two-dimensions  by  Ding  & Gidaspow 
(1990). 

4.2.  Fluidization  by  Lift  of  Spherical  Particles 

The  transport  of  particles  by  fluids  in  coal-water  slurries,  hydraulically  fractured  rocks 
in  oil-bearing  reservoirs,  bedload  transport  in  rivers  and  canals  and  their  overall  effect 
on  the  river  bed  erosion  etc.,  are  important  scientific  and  industrial  issues  in  particulate 
flows.  In  order  to  understand  fluidization/sedimentation  in  such  conduits,  (Choi  & Joseph 
2001)  performed  a DNS  study  of  fluidization  of  circular  cylinders  (300  particles)  arranged 
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Figure  3.  Temporal  evolution  of  particle  distribution  during  fluidization  by  a gas  jet.  Initially 
all  the  particles  are  uniformly  arranged  in  layers  at  the  bottom  of  the  rectangular  box.  Air 
is  injected  through  a rectangular  slot  at  the  bottom  wall.  Air  bubbles  are  trapped  within  the 
particles  and  the  growth  and  pattern  of  these  bubbles  are  in  agreement  with  experimental 
observations. 


Computational  domain,  63  x 12  x 12cm 
Gas  jet  velocity,  9m/s 
Fluid  density,  Ip/cm^ 

Particle  Density,  l.Olp/cm^ 

Numer  of  Parcels,  3780 
Initial  array  height,  4.75cm, 

Pressure  gradient,  20dj/ne/cm® 


Grid,  20  X 11  X 10 

Jet  diameter,  0.04m 

Fluid  viscosity,  Ipoise 

Diameter  of  particles,  0.95cm 

Particles  per  parcel,  1 

Initial  centerline  velocity,  360cm/s 


Table  3.  Parameter  description  for  the  simulation  of  fluidization  of  spherical  particles  in  a 

plane  Poiseuille  flow. 


at  the  bottom  of  a channel  in  plane  Poiseuille  flow.  They  observed  that  with  sufficient 
pressure  gradient  across  the  channel,  the  particles  initially  at  rest  in  the  lower  half  of  the 
channel  start  moving  and  roll  over  the  wall.  Particle  rotation  in  a shear  flow  generates  lift 
and  the  channel  is  fluidized  after  some  time.  We  attempt  to  capture  this  phenomenon  by 
our  Eulerian-Lagrangian  model.  The  flow  parameters  are  given  in  Table  3.  As  opposed  to 
Choi  & Joseph  (2001),  we  are  performing  three-dimensional  simulations  and  our  particles 
are  spheres.  As  shown  in  Fig.4  we  observe  that  the  effect  of  volume  displacement  due 
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Figure  4.  Temporal  evolution  of  particle  distribution  during  fluidization  by  lift  in  a plane 
Poiseuille  flow.  Also  shown  axe  contours  of  aocial  gas-phase  velocity.  Initially  all  particles  are 
arranged  at  the  bottom  of  the  channel. 


to  particles  is  to  set  up  Kelvin-Helmholtz  type  waves  between  the  fluid  and  particle 
layers..  A vertical  pressure  gradient  is  created  and  gives  vertical  velocity  to  the  particles 
and  the  channel  gets  fluidized.  We  also  did  several  test  cases,  with  higher  grid  resolution, 
increased  density  ratios  to  obtain  similar  results.  With  increased  particle  density,  the  two- 
way  coupling  between  the  particle  and  fluid  momentum  equation,  decelerates  the  fluid  in 
the  bottom  half  of  a channel  and  an  inflection  point  is  created  in  the  axial  velocity  profile. 
This  eventually  cause  lift  and  particle  dispersal.  In  liquid-fuel  combustor  applications, 
we  believe  that  similar  mechanism  can  be  observed  in  the  dense-spray  regime  near  the 
injector  and  will  allow  us  to  capture  the  important  features  of  primary  atomization. 


5.  Summary 

We  have  developed  a numerical  scheme  which  accounts  for  the  displacement  of  the  fluid 
by  particles  and  is  appflcable  to  dense  particulate  flows  and  spray  regimes  close  to  the 
nozzle  injector.  This  scheme  has  been  implemented  into  the  unstructured  LES  code.  We 
performed  various  numerical  simulations  of  fluidization  and  particulate  flows  to  validate 
our  scheme.  It  was  shown  through  simulation  of  channel  flow  with  layers  of  spherical 
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particles  at  the  bottom  wall  that  the  effect  of  particle  volume  fraction  is  to  displace  the 
fluid  in  a control  volume  and  create  vertical  pressure  gradient  through  two-way  coupling 
between  the  particles  and  fluid.  This  gives  Kelvin-Helmholtz  type  instability  waves  which 
result  in  dispersal  and  lift-off  of  particles.  This  liftoff  of  dense-particles  cannot  be  observed 
without  explicit  modeling  of  lift-forces  on  particles  by  using  the  standard  point-particle 
approximation  in  Lagrangian-Eulerian  framework.  Similar  K-H  waves  are  present  in  the 
dense  spray  regime  near  the  injector  nozzle  which  may  affect  the  overall  spray  dynamics 
and  breakup  in  practical  combustors.  We  plan  to  apply  this  model  coupled  with  secondary 
breakup  model  to  perform  simulations  of  spray  atomization.  We  also  plan  to  perform 
several  validation  studies  of  this  methodology  for  dense  particulate/granular  flows  to 
address  some  issues  related  to  collision  scheme,  efficient  and  consistent  computation  of 
particle  volume  fraction  on  multiple  domain-decomposition. 
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